1 /* Copyright 2002-2016 CS Systèmes d'Information
2 * Licensed to CS Systèmes d'Information (CS) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * CS licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17 package org.orekit.propagation.semianalytical.dsst.utilities.hansen;
18
19 import org.hipparchus.analysis.polynomials.PolynomialFunction;
20 import org.hipparchus.util.FastMath;
21
22 /**
23 * Hansen coefficients K(t,n,s) for t=0 and n < 0.
24 * <p>
25 *Implements Collins 4-242 or echivalently, Danielson 2.7.3-(6) for Hansen Coefficients and
26 * Collins 4-245 or Danielson 3.1-(7) for derivatives. The recursions are transformed into
27 * composition of linear transformations to obtain the associated polynomials
28 * for coefficients and their derivatives - see Petre's paper
29 *
30 * @author Petre Bazavan
31 * @author Lucian Barbulescu
32 */
33 public class HansenZonalLinear {
34
35 /** The number of coefficients that will be computed with a set of roots. */
36 private static final int SLICE = 10;
37
38 /**
39 * The first vector of polynomials associated to Hansen coefficients and
40 * derivatives.
41 */
42 private PolynomialFunction[][] mpvec;
43
44 /** The second vector of polynomials associated only to derivatives. */
45 private PolynomialFunction[][] mpvecDeriv;
46
47 /** The Hansen coefficients used as roots. */
48 private double[][] hansenRoot;
49
50 /** The derivatives of the Hansen coefficients used as roots. */
51 private double[][] hansenDerivRoot;
52
53 /** The minimum value for the order. */
54 private int Nmin;
55
56
57 /** The index of the initial condition, Petre's paper. */
58 private int N0;
59
60 /** The s coefficient. */
61 private int s;
62
63 /**
64 * The offset used to identify the polynomial that corresponds to a negative
65 * value of n in the internal array that starts at 0.
66 */
67 private int offset;
68
69 /** The number of slices needed to compute the coefficients. */
70 private int numSlices;
71
72 /** 2<sup>s</sup>. */
73 private double twots;
74
75 /** 2*s+1. */
76 private int twosp1;
77
78 /** 2*s. */
79 private int twos;
80
81 /** (2*s+1) / 2<sup>s</sup>. */
82 private double twosp1otwots;
83
84 /**
85 * Constructor.
86 *
87 * @param nMax the maximum (absolute) value of n coefficient
88 * @param s s coefficient
89 */
90 public HansenZonalLinear(final int nMax, final int s) {
91
92 //Initialize fields
93 this.offset = nMax + 1;
94 this.Nmin = -nMax - 1;
95 N0 = -(s + 2);
96 this.s = s;
97 this.twots = FastMath.pow(2., s);
98 this.twos = 2 * s;
99 this.twosp1 = this.twos + 1;
100 this.twosp1otwots = (double) this.twosp1 / this.twots;
101
102 // prepare structures for stored data
103 final int size = nMax - s - 1;
104 mpvec = new PolynomialFunction[size][];
105 mpvecDeriv = new PolynomialFunction[size][];
106
107 this.numSlices = FastMath.max((int) FastMath.ceil(((double) size) / SLICE), 1);
108 hansenRoot = new double[numSlices][2];
109 hansenDerivRoot = new double[numSlices][2];
110
111 // Prepare the data base of associated polynomials
112 generatePolynomials();
113
114 }
115
116 /**
117 * Compute polynomial coefficient a.
118 *
119 * <p>
120 * It is used to generate the coefficient for K₀<sup>-n, s</sup> when computing K₀<sup>-n-1, s</sup>
121 * and the coefficient for dK₀<sup>-n, s</sup> / de² when computing dK₀<sup>-n-1, s</sup> / de²
122 * </p>
123 *
124 * <p>
125 * See Danielson 2.7.3-(6) and Collins 4-242 and 4-245
126 * </p>
127 *
128 * @param mnm1 -n-1 value
129 * @return the polynomial
130 */
131 private PolynomialFunction a(final int mnm1) {
132 // from recurence Collins 4-242
133 final double d1 = (mnm1 + 2) * (2 * mnm1 + 5);
134 final double d2 = (mnm1 + 2 - s) * (mnm1 + 2 + s);
135 return new PolynomialFunction(new double[] {
136 0.0, 0.0, d1 / d2
137 });
138 }
139
140 /**
141 * Compute polynomial coefficient b.
142 *
143 * <p>
144 * It is used to generate the coefficient for K₀<sup>-n+1, s</sup> when computing K₀<sup>-n-1, s</sup>
145 * and the coefficient for dK₀<sup>-n+1, s</sup> / de² when computing dK₀<sup>-n-1, s</sup> / de²
146 * </p>
147 *
148 * <p>
149 * See Danielson 2.7.3-(6) and Collins 4-242 and 4-245
150 * </p>
151 *
152 * @param mnm1 -n-1 value
153 * @return the polynomial
154 */
155 private PolynomialFunction b(final int mnm1) {
156 // from recurence Collins 4-242
157 final double d1 = (mnm1 + 2) * (mnm1 + 3);
158 final double d2 = (mnm1 + 2 - s) * (mnm1 + 2 + s);
159 return new PolynomialFunction(new double[] {
160 0.0, 0.0, -d1 / d2
161 });
162 }
163
164 /**
165 * Generate the polynomials needed in the linear transformation.
166 *
167 * <p>
168 * See Petre's paper
169 * </p>
170 */
171 private void generatePolynomials() {
172
173 int sliceCounter = 0;
174 int index;
175
176 // Initialisation of matrix for linear transformmations
177 // The final configuration of these matrix are obtained by composition
178 // of linear transformations
179 PolynomialFunctionMatrix A = HansenUtilities.buildIdentityMatrix2();
180 PolynomialFunctionMatrix D = HansenUtilities.buildZeroMatrix2();
181 PolynomialFunctionMatrix E = HansenUtilities.buildIdentityMatrix2();
182
183 // generation of polynomials associated to Hansen coefficients and to
184 // their derivatives
185 final PolynomialFunctionMatrix a = HansenUtilities.buildZeroMatrix2();
186 a.setElem(0, 1, HansenUtilities.ONE);
187
188 //The B matrix is constant.
189 final PolynomialFunctionMatrix B = HansenUtilities.buildZeroMatrix2();
190 // from Collins 4-245 and Petre's paper
191 B.setElem(1, 1, new PolynomialFunction(new double[] {
192 2.0
193 }));
194
195 for (int i = N0 - 1; i > Nmin - 1; i--) {
196 index = i + offset;
197 // Matrix of the current linear transformation
198 // Petre's paper
199 a.setMatrixLine(1, new PolynomialFunction[] {
200 b(i), a(i)
201 });
202 // composition of linear transformations
203 // see Petre's paper
204 A = A.multiply(a);
205 // store the polynomials for Hansen coefficients
206 mpvec[index] = A.getMatrixLine(1);
207
208 D = D.multiply(a);
209 E = E.multiply(a);
210 D = D.add(E.multiply(B));
211
212 // store the polynomials for Hansen coefficients from the expressions
213 // of derivatives
214 mpvecDeriv[index] = D.getMatrixLine(1);
215
216 if (++sliceCounter % SLICE == 0) {
217 // Re-Initialisation of matrix for linear transformmations
218 // The final configuration of these matrix are obtained by composition
219 // of linear transformations
220 A = HansenUtilities.buildIdentityMatrix2();
221 D = HansenUtilities.buildZeroMatrix2();
222 E = HansenUtilities.buildIdentityMatrix2();
223 }
224
225 }
226 }
227
228 /**
229 * Compute the roots for the Hansen coefficients and their derivatives.
230 *
231 * @param chi 1 / sqrt(1 - e²)
232 */
233 public void computeInitValues(final double chi) {
234 // compute the values for n=s and n=s+1
235 // See Danielson 2.7.3-(6a,b)
236 hansenRoot[0][0] = 0;
237 hansenRoot[0][1] = FastMath.pow(chi, this.twosp1) / this.twots;
238 hansenDerivRoot[0][0] = 0;
239 hansenDerivRoot[0][1] = this.twosp1otwots * FastMath.pow(chi, this.twos);
240
241 final int st = -s - 1;
242 for (int i = 1; i < numSlices; i++) {
243 for (int j = 0; j < 2; j++) {
244 // Get the required polynomials
245 final PolynomialFunction[] mv = mpvec[st - (i * SLICE) - j + offset];
246 final PolynomialFunction[] sv = mpvecDeriv[st - (i * SLICE) - j + offset];
247
248 //Compute the root derivatives
249 hansenDerivRoot[i][j] = mv[1].value(chi) * hansenDerivRoot[i - 1][1] +
250 mv[0].value(chi) * hansenDerivRoot[i - 1][0] +
251 (sv[1].value(chi) * hansenRoot[i - 1][1] +
252 sv[0].value(chi) * hansenRoot[i - 1][0]
253 ) / chi;
254 hansenRoot[i][j] = mv[1].value(chi) * hansenRoot[i - 1][1] +
255 mv[0].value(chi) * hansenRoot[i - 1][0];
256
257 }
258
259 }
260 }
261
262 /**
263 * Get the K₀<sup>-n-1,s</sup> coefficient value.
264 *
265 * <p> The s value is given in the class constructor
266 *
267 * @param mnm1 (-n-1) coefficient
268 * @param chi The value of χ
269 * @return K₀<sup>-n-1,s</sup>
270 */
271 public double getValue(final int mnm1, final double chi) {
272
273 //Compute n
274 final int n = -mnm1 - 1;
275
276 //Compute the potential slice
277 int sliceNo = (n - s) / SLICE;
278 if (sliceNo < numSlices) {
279 //Compute the index within the slice
280 final int indexInSlice = (n - s) % SLICE;
281
282 //Check if a root must be returned
283 if (indexInSlice <= 1) {
284 return hansenRoot[sliceNo][indexInSlice];
285 }
286 } else {
287 //the value was a potential root for a slice, but that slice was not required
288 //Decrease the slice number
289 sliceNo--;
290 }
291
292 // Danielson 2.7.3-(6c)/Collins 4-242 and Petre's paper
293 final PolynomialFunction[] v = mpvec[mnm1 + offset];
294 double ret = v[1].value(chi) * hansenRoot[sliceNo][1];
295 if (hansenRoot[sliceNo][0] != 0) {
296 ret += v[0].value(chi) * hansenRoot[sliceNo][0];
297 }
298 return ret;
299 }
300
301 /**
302 * Get the dK₀<sup>-n-1,s</sup> / dΧ coefficient derivative.
303 *
304 * <p> The s value is given in the class constructor.
305 *
306 * @param mnm1 (-n-1) coefficient
307 * @param chi The value of χ
308 * @return dK₀<sup>-n-1,s</sup> / dΧ
309 */
310 public double getDerivative(final int mnm1, final double chi) {
311
312 //Compute n
313 final int n = -mnm1 - 1;
314
315 //Compute the potential slice
316 int sliceNo = (n - s) / SLICE;
317 if (sliceNo < numSlices) {
318 //Compute the index within the slice
319 final int indexInSlice = (n - s) % SLICE;
320
321 //Check if a root must be returned
322 if (indexInSlice <= 1) {
323 return hansenDerivRoot[sliceNo][indexInSlice];
324 }
325 } else {
326 //the value was a potential root for a slice, but that slice was not required
327 //Decrease the slice number
328 sliceNo--;
329 }
330
331 // Danielson 3.1-(7c) and Petre's paper
332 final PolynomialFunction[] v = mpvec[mnm1 + offset];
333 double ret = v[1].value(chi) * hansenDerivRoot[sliceNo][1];
334 if (hansenDerivRoot[sliceNo][0] != 0) {
335 ret += v[0].value(chi) * hansenDerivRoot[sliceNo][0];
336 }
337
338 // Danielson 2.7.3-(6b)
339 final PolynomialFunction[] v1 = mpvecDeriv[mnm1 + offset];
340 double hret = v1[1].value(chi) * hansenRoot[sliceNo][1];
341 if (hansenRoot[sliceNo][0] != 0) {
342 hret += v1[0].value(chi) * hansenRoot[sliceNo][0];
343 }
344 ret += hret / chi;
345
346 return ret;
347
348 }
349
350 }